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Numerical simulations of convection in a layer filled with ideal gas are presented. The control 
parameters are chosen such that there is a significant variation of density of the gas in going from 
the bottom to the top of the layer. The relations between the Rayleigh, Peclet and Nusselt numbers 
depend on the density stratification. It is proposed to use a data reduction which accounts for the 
variable density by introducing into the scaling laws an effective density. The relevant density is 
the geometric mean of the maximum and minimum densities in the layer. A good fit to the data 
is then obtained with power laws with the same exponent as for fluids in the Boussinesq limit. 
Two relations connect the top and bottom boundary layers: The kinetic energy densities computed 
from free fall velocities are equal at the top and bottom, and the products of free fall velocities and 
maximum horizontal velocities are equal for both boundaries. 
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I. INTRODUCTION 

There is by now a large body of literature on turbulent 
Rayleigh-Benard convection [1] . Most of this work strives 
to stay in a limit described by an approximation devel- 
oped by Oberbeck and Boussinesq, commonly known as 
the Boussinesq approximation [2|. Within this approx- 
imation, the material constants are uniform across the 
layer and the temperature difference between top and 
bottom boundaries is small compared with the absolute 
temperature at any point in the layer. This idealized sys- 
tem has served as a paradigm for more complicated con- 
vection problems, as for example atmospheric convection. 
A major difference between convection in the Boussinesq 
limit and convection in the atmosphere is that the latter 
occurs in a layer in which the gas density varies signifi- 
cantly, which implies concomitant variations of viscosity 
and thermal diffusivity. 

Experiments on convection beyond the Boussinesq ap- 
proximation have mostly focused on effects caused by 
the temperature dependence of the material properties 
0, 01 J ^ opposed to the effects of compressibility. Exper- 
iments in low temperature gases near their critical point 
are to some extent an exception because the parameters 
in these experiments can be adjusted such that the adi- 
abatic temperature gradient in the gas delays the onset 
of convection and shifts the critical Rayleigh number by 
a detectable amount |5|-i7||. However, these experiments 
were restricted to convection near the onset. Experimen- 
tal studies aimed at turbulent Boussinesq convection in 
low temperature gases also have to correct for the effect 
of the adiabatic temperature gradient [1, Q . 

The present paper investigates through numerical sim- 
ulation convecting ideal gas in a layer with significant 
density variation from top to bottom. No slip top and 
bottom boundaries are employed, so that the results are 
in principle amenable to verification by laboratory ex- 
periments if one finds a way of realizing similar density 
gradients in turbulent convection experimentally. The 



parameters controlling the density variation and the adi- 
abatic temperature gradient are chosen to scatter around 
the values of the terrestrial troposphere. The troposphere 
is the bottom layer of the atmosphere of approximately 
10 km thickness which is well mixed by convection and 
bounded from above by the stably stratified stratosphere 

A long standing question in turbulent Rayleigh-Benard 
convection has been how the heat transport across the 
layer depends on the control parameters. One goal of the 
present paper will be to find out whether known results 
on convection in an incompressible medium can be ex- 
tended to an ideal gas. From a fundamental point of view, 
the scale height of the density profile introduces a new 
length scale into the problem in addition to the height of 
the layer. Previous studies [ll| suggest that convective 
motion extends through multiple scale heights, so that 
we cannot expect that the layer height will simply drop 
out of the list of the relevant parameters in order to be 
replaced by the scale height. 

Another issue peculiar to non-Boussinesq convection is 
the asymmetry between the boundary layers next to the 
top and bottom boundaries. For example, the heat con- 
ductivities near the warm bottom and cold top bound- 
aries are different, but the heat fluxes through both 
boundaries are identical in a statistically stationary state, 
so that the temperature gradients within the top and 
bottom boundary layers must be different. Because of 
this asymmetry, the temperature in the center of a con- 
vection cell does not need to be equal to the arithmetic 
mean of top and bottom boundary temperatures. The 
deviation of the true center temperature from this arith- 
metic mean has been used in experiments as an indicator 
of non-Boussinesq effects [3, |3| • 

A typical situation in astro- and geophysics is that only 
part of a convective layer is accessible to observation, 
at least to accurate observation. For example, the top 
of planetary or stellar atmospheres and the bottom of 
Earth's troposphere are better known than the rest of 
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the convective layers. It is in these cases important to 
know what can be inferred about the convective layer 
from the observation of some part of it. Translated to 
the idealized system simulated here, the question arises 
as to what can be deduced about the whole convective 
layer or a boundary layer from knowledge of the opposite 
boundary layer. 

The next section will present the mathematical model 
and the numerical method used to solve it. The numeri- 
cal results are analyzed in the third section, in which one 
subsection deals with the scaling of heat transfer and 
kinetic energy with the control parameters, whereas an- 
other subsection is concerned with the relationship be- 
tween the boundary layers. 



thermodynamics that R = Cp — Cv- The strain rate 
tensor is given by = ^{djVi + diVj). 

From here on, we will use nondimensional variables. 
All lengths are expressed in multiples of d, and the scales 
of time and density are chosen as (P / Ko and po, respec- 
tively. The difference between the temperature of the 
gas and the top temperature, T, is scaled with AT. Us- 
ing the same symbols for the non-dimensional variables 
space, time, density, velocity and temperature difference 
with the top boundary as in the dimensional equations 
([T]|4]), one obtains the system 

dtp + V- (pv) - (5) 



II. MATHEMATICAL MODEL AND 
NUMERICAL METHODS 

Consider a plane layer of height d bounded by two 
planes perpendicular to the 2;— axis. Gravity g is con- 
stant and pointing along the negative z— axis, g = ~gz, 
where the hat denotes a unit vector. The ideal gas is 
characterized by constant heat capacities at fixed volume 
and pressure, Cy and Cp, and constant dynamic viscos- 
ity /i and heat conductivity k. This implies that the 
density dependences of kinematic viscosity 1/ and ther- 
mal diffusivity k are given hy fi — pv and k — upCp, 
in which p is the density. Let us assume that top and 
bottom boundaries are no slip and have prescribed tem- 
peratures. Parameters evaluated at the top boundary in 
the initial state will be denoted by an index o for "outer" . 
The gas at the top of the layer thus has temperature Tq. 
In the state specified by the initial conditions (see Eq. 
()13p below), it also has kinematic viscosity thermal 
diffusivity Kq and density po ■ The temperature difference 
across the layer is AT. 

The system of equations governing density, tempera- 
ture T -\- To, pressure p and velocity v reads with the 
usual summation convention over repeated indices: 



dfV + {v ■ \7)v = 



-^V[(T+|^)p]i^PrRa (6) 



zPt Ra 



AT 



3 p 



dtT + v-VT = 



lv^T~i^-l)iT+^)V-v (7) 



1 



.1 



d 1 



together with the boundary conditions 



Ho Ra 



r(z = l)=0 , T(z = 0) = l , v{z = l)=v{z = 0)=0. 

(8) 

The equation of state has been used to eliminate pressure. 

Seven parameters control the system. The Rayleigh 
number Ra is the usual Rayleigh number evaluated at 
the top boundary (remember that the thermal expansion 
coefficient of an ideal gas is its inverse temperature) : 



Ra 



gd^AT 

ToK 0^0 



(9) 



dtp + V- {pv) = 



(1) 



The Prandtl number Pr is independent of space in the 
present model and is set to 0.7 in all calculations: 



p[dtv + {v-V)v] = -Vp + pg + p[V^v + -V{V-v)] (2) 



dtT+vVT 



p 



Cl^^'^ pCy 



V-v- 



2p 



p = pR{T + To) 



(3) 



(4) 



The gas constant R in the equation of state is given 
by i? = Ru/m, with the molar mass m and the universal 
gas constant i?„ = 8.314J mol~^ . It follows from 



Pr=^ 

K 



0.7. 



(10) 



The adiabatic exponent 7 is set in this paper to its value 
for a monoatomic gas: 



7 = Cp/Cv = -. 



(11) 



The density stratification is specified by d/Ho, where 
Ho is the adiabatic scale height at the top boundary. 
Ho — ^RTo/g. The meaning of the fifth parameter, 
AT /To, is obvious from the definitions above. An alter- 
native parameter, redundant after the choices made so 
far, is the ratio of the adiabatic temperature difference 
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between top and bottom, ATad, 
ature difference, AT: 



and the actual temper- 



(a) 



AT, 



ad 



AT 



(7-1) 



ATH, 



gd 



CpAT 



(12) 



This ratio needs to be less than 1 for any convection to 
occur. 

The sixth "parameter" is the initial temperature and 
density distribution. The initial conditions appear as a 
control parameter because they specify for instance the 
total mass in the layer. They also determine po and hence 
Vo and Ko, a quantity used to make the governing equa- 
tions nondimensional. All simulations are started from 
zero velocity and the conductive profile T = 1 — z. The 
density is then determined from ([B]): 



1-71 



P = Po 



1 



(13) 



To /AT , 



The geometry, quantified through the aspect ratio of 
the computational volume, is the seventh parameter. 
Periodic boundary conditions are imposed in x— and 
y— directions with periodicity lengths Ix and ly. All com- 
putations have been made for Ix = ly = 2(i. No aspect 
ratio dependence has been investigated since the main 
interest of the present work was to determine the effects 
of density variations. 

Even though it was chosen to keep several parame- 
ters fixed, there still remains a vast parameter space to 
explore since Ra, d/Ho and AT /To have to be varied. 
The computations below are roughly guided by the ter- 
restrial troposphere [13], for which AT/Tq k, 0.35 and 
d/Ho ~ 1.9. The density of air varies by a factor be- 
tween 3 and 4 within the troposphere. Note that 7 = 1.4 
is the appropriate adiabatic exponent for air. 

The well known Boussinesq approximation is recov- 
ered from equations ([5][7]) in the limit d/Ho — >■ and 
AT /To — >■ 0. In this hmit, the sound speed goes to 
infinity. A purely explicit time step will be used in the 
numerical method below, so that simulations close to the 
Boussinesq limit become impractical. For this reason, 
and since the Boussinesq limit is an important reference 
case, a second system of equations was also implemented 
numerically: 



(14) 



dtv + {v ■ W)v = -c^Wp + Pr Ra Oz + PrW^v (15) 



(16) 



where 9 is the deviation from the conductive profile, 
T = 1 — z + 9. The sound speed c appears as an inde- 
pendent variable. For small Mach numbers, the density 
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FIG. 1: (Color online) Variations as a function of height z of 
horizontal averages of temperature T (top panel), density p 
(second panel from the top), vertical velocity Vz (third panel) 
and horizontal velocity (bottom panel). The overbars signal 
averages over time and horizontal planes. The different traces 
are for — = 3 and the Rayleigh numbers 2 x 10^ (solid 
red line), 2 x 10^ (long dashed green line), 2 x 10^ (dot dashed 
blue line), 2 x 10^ (short dashed black line). 
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fluctuations are small and the equation of continuity can 
be linearized to yield (|14l) . The Boussinesq equations 
are recovered in the limit c — > oo. In the simulations 
mentioned below using (|14m6p . the sound speed was ad- 
justed such that the Mach number always stayed below 
0.1. Note that simulations of weakly compressible con- 
vection as an approximation to Boussinesq convection 
have been undertaken before. For instance, simulations 
using the lattice Boltzmann method [l^ implicitly do so. 

Systems ([MZl) and (|14m6p have been simulated with a 
finite difference method implemented on graphic process- 
ing units using C for CUDA. The numerical method used 
centered finite differences of second order on a collocated 
grid except for the advection terms which used an upwind 
biased third order scheme. The time step was a third or- 
der Runge-Kutta method. The standard resolution was 
256'^. Lower resolution was sufficient at the smallest sim- 
ulated Ra. The validation of the code is described in the 
appendix. 



III. RESULTS 

A. Overview 

A summary of the simulations is given in tableHl Apart 
from the control parameters, it lists the Nusselt number, 
defined as 



Nu 



dT 

dz 



(17) 



where the overbar denotes average over time and either 
top or bottom boundary. The kinetic energy density i?kin 
is given by 



I f 1 



-pv^dV 



V J 2' 

whereas the Peclet number Pe is computed from 
Pc = 




(18) 



(19) 



which is aptly called Peclet number because velocities 
are computed in units of Ko/d. The average tempera- 
ture deviation from the conductive profile at the center 
of the layer, O^m is also listed in table U together with 
the average density in the midplane, pm- 

Fig. [1] shows vertical profiles of temperature, density, 
vertical velocity, and horizontal velocity (u^ -t- VyY^"^ for 
different Ra. Contrary to the Boussinesq case, these pro- 
files are not symmetric about the midplane. The relation 
between the top and bottom regions of the layer will be 
discussed in section IIIICl As Ra increases, an increas- 
ingly large interval develops in which the temperature 
gradient is approximately equal to the adiabatic gradi- 
ent. The maximum of the vertical velocity is found below 
the midplane. Two local maxima show up in the profiles 



Nu 




FIG. 2: (Color online) Nusselt number Nu as a function 
of Rayleigh number Ra within the Boussinesq approxima- 
tion (black stars), for ATad/AT = 1/15 (blue symbols) 
and at/To = 0.1 (plus) and 1 (x), for ATad/AT = 2/3 
(green symbols) and AT/To = 0.1 (empty squares), 0.3 
(full squares), 1 (empty triangle up), 3 (full triangle up), 
10 (empty triangle down), 30 (full triangle down) and 100 
(empty diamonds) and for ATad/AT = 4/5 (red symbols) 
and at/To = 0.1 (empty circles), 1 (full circles) and 10 (half 
filled circles). Data for the Boussinesq case do not appear in 
the subsequent figures. 



of horizontal velocity. The larger velocities are found 
near the bottom boundary. When both d/Ho and Ra 
are large, there is no local maximum of horizontal veloc- 
ity near the top boundary and the horizontal velocity de- 
creases monotonically with height. These cases result in 
blanks in the last column of table HI Boundary layers also 
exist in the density profiles. The average density takes 
its maximum and minimum values, Pmax and Pmin, near 
the bottom and top boundary, respectively. Both pmax 
and pmin are given in table HI too. The ratio Pmax/Pmin 
is generally less than is suggested by the values of d/Ho 
and the initial conditions because the statistically 
stationary, turbulent and well mixed state is nearly adi- 
abatic, not conductive. A rough estimate of Pmax/Pmin 
can be obtained from assuming that the adiabatic state 
extends throughout the layer, implying that the density 
takes its extremal values exactly on the boundaries. This 
leads to 



Pn 



T 



AT, 



ad 



AT, 



ad 



1/(7-1) 



(20) 



with T„ 



(Om + 5) AT + To, which is compatible with 
the numerical results. 



B. Global Quantities 

The most immediate task is of course to find predic- 
tions for the Nusselt number. A straightforward plot of 
Nu vs. Ra (fig. [5]) shows that one does not obtain simple 
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(a) 10' 



FIG. 3: (Color online) Nu* as a function of Ra» with the same 
symbols as in fig. [2] The solid line indicates the power law 



Nu. = Ra'iV8-5. 



power laws for large ATad/AT and d/Ho- A finite adi- 
abatic temperature gradient modifies the onset of con- 
vection. If one aims for a data reduction which collapses 
the different curves in fig. [21 one can account for the 
adiabatic temperature difference by defining a corrected 
Rayleigh number Ra, by 



Ra* 



To Kq^O 



= Ra 1 - (7 - i; 



d To 



HoAT^ 

(21) 

A similar correction seems in order for Nu. Since the 
adiabatic temperature gradient needs to be established 
before any convection can start, it is natural to subtract 
the heat conducted down the adiabat from both the ac- 
tual heat transport and the conductive heat transport 
used for the normalization of the heat transport [isj : 



Nu* = 



Nu - ATad/AT 



(22) 



f - ATad/AT ■ 

It is seen from fig. [3] that for Ra* larger than roughly 
10^, one finds approximately (Nu» — 1) cx Ra°'^, but the 
prefactors depend on the other control parameters in a 
non-trivial way. This is compatible with an argument ex- 
posed in Ref. [ij, which states that Ra(Nu— 1) = /(Ra,) 
with / an a priori unknown function. The argument has 
to assume negligible viscous heating and small density 
variations, so that it is not expected to hold through- 
out the parameter range investigated here. Nonetheless, 
a best fit to the data yields Ra(Nu - 1) = Ra^'^/S.S 
which can be rearranged into the fitting function in fig. 



1 



Ra' 



0.3 



m (Nu, 

It can also be useful to relate Nu to the kinetic en- 
ergy or the Peclet number. It was noted in Ref. [l^ that 
(Nu — I) oc Pe^^'^ in Boussinesq convection in computa- 
tional volumes of large aspect ratio, which is equivalent to 

1 /3 

(Nu-1) OC E^r in that case. In the present simulations, 
the relation between i?kin and Pe is already non-trivial 
(see fig. S]), because there is a factor representing an ef- 
fective density between the two quantities. It turns out 
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FIG. 4: (Color online) The kinetic energy density -Ekin as 
a function of Pe (top panel) and Pe(pmaxPmin)^'^^ (bottom 
panel) with the same symbols as in fig. [2l The solid 
lines are given by i?kin = ^Pe^ (top panel) and i?kin = 



^VPmaxPmin Pe^ (bottom paucl) . 



that the geometric mean of Pmax and pmin is a suitable 
effective density to the extent that in fig. SI all points 
for Ra,(pinaxPmin)^^^ > 100 deviate by less than 30 % in 
iJkin from 



-^kin — 2 's/ PmaxPmin P^ 



(23) 



This geometric mean becomes again important when 
looking for a relation between Nu, — 1 and i^kin- A good 
fit to the data is obtained from 



Nu, 



1 — y {E]iin\/ PmaxPn 



Nl/3 



(24) 



(see fig. [S]) which reduces of course to the previously 
known scaling [H, [11] for pmax = Pmin and ATad/AT = 
0. 

Having established the relevance of the product 
PmaxPmin, it is tempting to introduce it into fits of Nu, 
vs. Ra,. A reasonable fit is shown in fig. [6] to be 

Nu, -0.22(Ra,(p„,axPmin)'/^) , (25) 

which is an improvement compared with fig. |3] especially 

for Ra,(pmax/Omin)^^'' > 10''. 
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FIG. 5: (Color online) Nu, — 1 as a function of 
i5kin(pmaxPmiii)^'''^, Compensated for the power law in Eq. 
(|24|l . with the same symbols as in fig. [2] 



C. Boundary Layers 

Previous studies have quantified the asymmetry be- 
tween top and bottom in non-Boussinesq convection with 
the help of the midplane temperature [1, 01 • In most of 
the present simulations, the midplane temperature de- 
viates by less than 0.05 AT from its value in the con- 
ductive state (see table Such a small deviation is 
difficult to determine accurately and requires long time 
integrations, so that this section will not consider 6m any 
further, apart from noting that 9m is negative for small 
ATad/AT (in agreement with Ref. H) but becomes posi- 
tive for ATad/AT large enough. 

A relation between temperature boundary layers de- 
duced from experimental data by Wu and Libchaber Q 
is based on a temperature scale computed from quantities 
local to each boundary layer. In many of the simulations 
presented here, the boundary layers are still quite thick 
and there is significant variation of for example thermal 
diffusivity across them, so that the results of Wu and 
Libchaber cannot be tested in a meaningful way. 

In the following, overbars denote averages over time 
and horizontal planes, and the indices b and t indicate 
top and bottom boundaries. For example, pt is the aver- 
age density at the top boundary. It is in general different 
from po, which is the density at the top boundary in the 
initial conductive state given by Eq. (fT3|) . 

This subsection will present two relations between the 
top and bottom regions of the convection layer involving 
the free fall velocity. In order to compute a free fall 
velocity, we define Tad,t and Tad.b as the (dimensional) 
temperatures the gas would have at the top and bottom 
boundaries if the adiabatic temperature profile extended 
throughout the entire layer: 



Tad,t — Tm — — ATad 



-^ad,b ^rr. 



-ATad 

2 



(26) 



with T„ 



'i)AT-l-To. Under the same assumption. 
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FIG. 6: (Color online) Nu, as a function of Ra*(pniaxPmin)^''*, 
on a double logaritmic scale in the top panel and compensated 
for the power law in Eq. (|25[) in the bottom panel, with the 
same symbols as in fig. [2l 



the densities at the two boundaries, Pad.t and Pad.b, are 
given by 



Pad.t 



Tad A 

T 



1/(7-1) /rj. X 1/(7-1) 



Pad,b 



/ Tad,b 



(27) 

Consider now a parcel of gas near the top boundary. It 
has on average the density pt- The density difference 
with the adiabatic profile, pt — Pad.t, accelerates the par- 
cel through the volume. The free fall velocity is esti- 
mated from a balance between the advection and buoy- 
ancy terms, which reads in the non-dimensional variables 
used here \p{v ■ V)v\ ^ Apgd{d/ Ko)^ , where Ap is the 
density difference of the moving parcel with the adia- 
batically stratified background. The pressure variation 
experienced by the falling parcel compresses the parcel 
by the same factor as the surrounding gas (assuming the 
parcel does not exchange heat with its surroundings) , so 
that Ap/p remains constant during the entire journey 
through the adiabatically stratified layer. It follows that 
Ap/p keeps its initial value of (pt — Pad,t)/Pt, and that 
the square of the non-dimensional free fall velocity of 
the parcel arriving at the bottom (which is expressed in 
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FIG. 7: (Color online) -Ekin as a function of Pad,b«ff,b with the 
same symbols as in fig. [21 The solid line indicates the power 
law £ki„ = 0.12 (pad,b4,b)° ''^ 
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FIG. 8: (Color onhne) (pad,t'ui,t)/(Pad,bt'ff.b) as a function of 
Pad.b'i^ff b with the same symbols as in fig. [21 



units of {Ko/(£f) is ^* p^"^'^ 9'^ (r") ■ ^ similar expres- 
sion is derived if we start the argument from the bottom 
boundary, so that we obtain two velocities, wg.t and fff.b 
according to the formula 



Pt — Pad.t r, r, "^o 

v., = 1 — ^^Pr 



1/2 



1/2 



V Pb 



AT J 



(28) 



Fig. |7| verifies that one obtains with Eq. a veloc- 
ity representative of the convective velocity. The figure 
shows -Ekin as a function of the energy density computed 
from the bottom free fall velocity, pad.bf^j,. When this 
velocity is small, the Reynolds number of the flow is too 
small for friction to be negligible and the free fall velocity 
is a poor estimate of the true velocity. For large veloci- 
ties, there is a unique relation between Ekin a-nd Pad.b^^i b 
independent of any other parameters. 



FIG. 9: (Color onhne) (wff ,tt^h,max,t)/(t^ff,b«h,max,b) as a func- 
tion of Wfi,bt'h,max.b with the Same symbols as in fig. [21 
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FIG. 10: (Color online) Wfi,b/«h,max,t as a function of Dh,max,t 
with the same symbols as in fig. [21 



It is now expected that one obtains the same graph 
if one uses the kinetic energy density computed at the 
top boundary, or equivalently, that Pad,tVg ^ = Pa.d.bVg b- 
Fig. [8] demonstrates that this is the case to within ±20% 
for sufficiently large velocities and over three decades in 
Peid,hVg b- The equality of the two energy densities is the 
first important connection between the top and bottom 
boundaries. 

The second relation, shown in fig. |S1 involves the local 
maxima of horizontal velocity near the top and bottom 
boundaries visible in the bottom panel of fig. |T] These 
maximum velocities, Uh,max,t and 'i'h,max,b, are listed in 
table m According to fig. ^ they obey for sufficiently 
large velocities: 



fff,tl'h.max,t = t^ff,bWh,max.b- 



(29) 



^'h,inax,b IS a free fall velocity computed from quantities 
evaluated at the the bottom boundary, but it estimates 
the velocity of plumes arriving at the top boundary. The 
typical velocity of the flow out of the arriving plumes is 
the maximum average horizontal velocity, so that we ex- 
pect i;h,max,t oc vg,h, and similarly Wh,max,b oc V{f,t. The 
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prefactors depend on the size of the incoming plumes 
relative to the thickness of the layer in which the out- 
flow occurs. Fig. [TOl shows that "Dh^max.t oc vg^h approx- 
imately holds and that the prefactor indeed depends on 
d/Ho and AT/Tq. The two proportionalities combine 
to Wff,t'i'h,max,t c< f ff,bi^h,inax,b- It is remarkable that this 
combined relation is obeyed more accurately than the 
two separate proportionalities, and that the proportion- 
ality factor in the combined relation approaches 1 at high 
Rayleigh numbers as shown in fig. |9l 



IV. CONCLUSION 

There are many ways how to depart from the Boussi- 
nesq approximation 0, 0, @, 0] and it is not known at 
present whether there is anything universal about con- 
vection in general. The Prandtl number is a constant in 
a real gas and the Nusselt number depends only on the 
Rayleigh number in the Boussinesq limit. Away from 
that limit, several more control parameters determine 
the Nusselt number and it is a challenge to find a suit- 
able data reduction such that the Nu{Ra) dependences 
found for various density stratifications collapse on a sin- 
gle master curve, which then must contain the Boussinesq 
limit. It has been shown in the present paper that a good, 
if imperfect, data collapse is obtained if one introduces 
into the scaling laws an effective density given by the ge- 
ometric mean of the maximum and minimum densities 
in the convecting layer. There is no theory underpinning 
the relevance of the geometric mean. It seems likely that 
the geometric mean will be of no help in some other cases 
of non-Boussinesq convection, for example in liquids. 

On the other hand, the importance of free fall velocities 
has been recognized since long, and it has been shown 
here that two different estimates of free fall velocities, 
based on quantities pertaining either to the top or the 
bottom boundary, are related at high Rayleigh numbers 
by the requirement that the kinetic energy densities com- 
puted from the two velocities must be equal. The free fall 
velocities are also connected to the total kinetic energy 
of the flow and the maxima of the horizontal velocity 
profile. It will be worthwhile to check these relations in 
other forms of convection. 



Appendix 

This appendix describes a few tests that have been 
performed in order to validate the numerical code. Direct 
validation of the code is problematic because published 
simulations of compressible convection either focus on 
flow structures [lllllTj and are not useful as a benchmark, 
or are too close to the Boussinesq limit to offer a stringent 
test [l^. A different route had therefore to be taken. 

First, the code simulating (|14j|16p could be validat;ed 
against a completely independent spectral method 19]. 
This already verifled most terms occurring in (l5][71). 



For example, the spectral code calculates for a plane 
layer with no slip boundaries in z, periodic bound- 
ary conditions in x and y with Ix = ly = '2d, and 
Pr = 0.7 and Ra = 1.6 x 10^ that Skin = 360 and 
Nu = 3.01. If started from the initial conditions v = 0, 
9 — sin(7rz) cos(27ry) cos(27rx) a,t t — 0, the system pH - 
[TB)) needs to be integrated beyond i = 10 to find the final 
state. The new code with appropriately adapted bound- 
ary conditions yields i^kin = 351 and Nu = 2.96 for 32 
points along z and 64^ grid points in the s,?/— plane. 
At twice that resolution in each coordinate, the result is 
iiJkin = 357 and Nu = 2.99. The parameter in Eq. 
((T5)) was set to 10^ so that the Mach number is below 
8.5 X 10-2. 

In a next step, the code implementing ([S][7]) was used 
to simulate the propagation of sound waves. For this 
purpose, the term —zPr Ra-^ was removed from Eq. 
([6|). The remaining equations, if linearized and neglecting 
dissipation, can be manipulated into 



AT d 



^Pr Ra V(V • v). 



This equation allows some simple analytical solutions. 
For example, there are the eigenmodes depending only on 
z. The first of those eigenmodes for boundary conditions 
imposing zero heat flux at the top and bottom boundaries 
is of the form Vz oc sin(7rz) cos(a;t). This standing wave 
has a period of t = 2tt/uj = 2{^^Pt Ra)-^^^. For 
4f = 0.1, = 0.01, Pr = 0.7 and Ra = 2 x 10^, this 
predicts r = 5.345 x 10^"*. With a resolution of 64 grid 
points along z, the numerical result is r = 5.343 x 10^*. 
One can similarly simulate sound waves propagating in 
different directions. This is a test for all terms involving 
V • t>. The fact that simulations of convection at high 
Ra yield a temperature gradient in the bulk close to the 
adiabatic gradient may also be regarded as a test of the 
terms involving V • v. 

The dissipation rate of the standing wave of the previ- 
ous paragraph can be used to test the dissipative term. 
A more complete test is provided by the energy budget 
which also tests the viscous heating term in the tempera- 
ture equation. If one takes the scalar product of Eq. ([5]) 
with V and integrates Eqs. ^ and ([7]), multiplied by p, 
over all space, one deduces from Eqs. that the time 

derivative of kinetic plus internal energy is given by 



d_ 

di 

with 



1 , 



G 

V2 



Ra Pr 

d J J 



To 



-{T- 



at' 



dV^G+Vi+V2 



-Pr Ra^ / pVzdV, 



Pr I v,[W\, + ^d,{V ■ v)]dV, 
Pr J 2[e^j-^V -vS^jfdV. 
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G is the work done by gravity, Vi the dissipated ki- 
netic energy and V2 the heat generated through viscous 
dissipation. If we denote time averages by angular brack- 
ets, we must find in a statistically stationary state that 
(G) = and {Vi) + {V2) = 0. Since Vi and V2 have 
different forms and must be programmed differently, the 
energy budget provides a good test for their correctness. 
For the case ~ 'n~ ~ ^ = 6 x 10^ included in 

tableU the simulations yield | (G) /{Vi)\^ 3.3 x 10"^ ^nd 
I (Vi V2)/(Vi) I = 7.6% at a resolution of 64^ grid points, 
and \{G)l{Vi)\ = 4.2 x lO'^ and K^i + V2) / {Vi)\ = 4% 



at a resolution of 128^ grid points. The formulae for Vi 
and V2 taken together contain integrals of 18 derivatives, 
6 of which are squared, so that the typical error on each 
derivative is about 0.1%. The volume integrals have been 
computed by adding the integrands at each grid point, 
multiplied by the volume of the cell surrounding each 
grid point. This method of integration is of first order 
for general integrands [20], which explains why the er- 
ror in \ {Vi + V2)/(Vi)| is only halved when doubling the 
resolution. 
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TABLE I: Summary of results. Listed are the control param- 
eters -i-, and Ra together with Nu, Skin, Pe, Om multi- 
plied by 100, pm, pmax, pmin, «h,max,b, and «li,max,t (seC toxt 

for definitions). The table consists of three sections (corre- 
sponding to the color code of the figures in the online version) 
with diflferent ATad/AT, which is 4/5, 2/3 or 1/15 in going 
from the top to the end of the table. An entry is missing for 
Vh,max,t if the profile of horizontal velocity has no maximum 
near the top boundary. 
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TABLE II: Table [I continued. 



